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Using molecular dynamics simulation we examine changeovers among crystal, glass, and liquid at 
high density in a two dimensional binary mixture. We change the ratio between the diameters of 
the two components and the temperature. The transitions from crystal to glass or liquid occur with 
proliferation of defects. We visualize the defects in terms of a disorder variable Dj (t) representing a 
deviation from the hexagonal order for particle j. The defect structures are heterogeneous and are 
particularly extended in polycrystal states. They look similar at the crystal-glass crossover and at 
the melting. Taking the average of Dj(t) over the particles, we define a disorder parameter D(t), 
which conveniently measures the degree of overall disorder. Its relaxation after quenching becomes 
slow at low temperature in the presence of size dispersity. Its steady state average is small in crystal 
and large in glass and liquid. 

PACS numbers: 61.43.-j, 61.72.-y, 61.70.Pf 



I. INTRODUCTION 

The phase behavior of binary particle systems is much 
more complicated than that of one component systems, 
where the temperature T, the average number density 
n = N/V, and the composition are natural control pa- 
rameters. At high densities, it is known to be profoundly 
influenced also by the size ratio <j\ / <ii between the diam- 
eters of the two components, g\ and oi [lj, |2j, y, |5( . If 
cti/(T2 is close to unity at large n, the system becomes a 
crystal at low T or a liquid at high T. If <7\/<J2 consider- 
ably deviates from unity, glass states are realized at large 
n and at low T. In glass states, the particle motions are 
nearly frozen and the structural relaxation time grows, 
but the particle configurations are random yielding the 
structure factors similar to those in liquid. 

Recently, the liquid-glass transition has been studied 
in a large number of molecular dynamics simulations on 
model binary mixtures both in two and three dimensions 
[3. 17L 13. EL llflj . In these simulations, the temperature T 
has mostly been the control parameter at fixed average 
density and composition. Some authors have applied a 
shear rate or a stress to glassy systems as a new control 
parameter |^, . The size ratio o\ / a 2 has been cho- 
sen at particular values to realize fully frustrated parti- 
cle configurations and to avoid crystallization and phase 
separation. However, for weaker size dispersity, the de- 
gree of disorder should become smaller. Polycrystals will 
be realized at some stage and a crystal with a small num- 
ber of point defects will be reached eventually. On this 
crossover we are not aware of any systematic study and 
have no clear picture. 

In this paper, we first aim to visualize the disorder 
brought about by the size dispersity in two dimensions 
(2D). To this end we will introduce a disorder variable 
Dj (t) > representing a deviation of the hexagonal crys- 
tal order around each particle j. Snapshots of Dj(t) re- 
alized by each simulation run will exhibit patterns indi- 
cating the nature of the defect structure. We shall ob- 



serve point defects in crystal, grain boundaries in poly- 
crystal, and amorphous disorder in glass. The average 
D(t) = J^jLi Dj(t)/N over the particles is a single "dis- 
order parameter" characterizing the degree of overall dis- 
order. 

Halperin and Nelson found that defects play a key 
role in 2D melting in one component systems, predict- 
ing continuous transitions with an intermediate " hexatic" 
phase between crystal and liquid. They introduced a 
sixfold orientation order variable, written as Xj m this 
paper. The correlation function g§ (r) of the thermal 
fluctuations of \j has been used to characterize the 2D 
defect-mediated mel ting th eoretically 0, 0, ^| and ex- 
perimentally [13|)[lJ)[lSll6j- Our disorder variable Dj(t) 
will be constructed from their \ji so we W1 h visualize 
the defect patterns exhibited by Dj(t) also at the melt- 
ing. The problem becomes much more complex for binary 
mixtures, where the crystal-liquid transition occurs with 
changing T or n at weak size dispersity 0, and the 
glass-liquid transition occurs at stronger size dispersity 
@ S IE HI Ejl ■ We should understand the defect struc- 
ture by changing <J\jai and T (and/or n) both at the 
crystal-glass and crystal- liquid transitions |17| . 

In Sec II, we will introduce the quantities mentioned 
above and present our numerical results at fixed den- 
sity and composition, where the defects involved in the 
crystal-glass and crystal-liquid transitions will be visual- 
ized. We will also calculate the overall disorder param- 
eter D(t) in transient states and in steady states as a 
function of <J\j<Ji and T . In Sec III, we will summarize 
our results and give some remarks. 



II. NUMERICAL RESULTS 
A. Method 

We used a 2D model binary mixture interacting via 
a truncated Lenard- Jones (LJ) potential v a p{r\ where 
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a, j3 — 1,2 represent the particle species. If the distance 
r between two particles is larger than a cut-off r cut , we set 
v a p(r) = 0. If r < r cut , it is given by the Lennard- Jones 
potential, 



v a p(r) = 4e 



(^f -(- 



a 



a/3, 



(1) 



which is characterized by the energy e and the soft-core 
diameter a a p — (a a + ap)/2 with <j\ and a-i representing 
the (soft-core) diameters of the two components. The 
constant C a p ensures v a p{r) — > as r — > r cut , so the 
potential is continuous at the cut-off distance. We set 
''cut = 3.2f7i for any a and j3 [Isj . The particle numbers 
of the two species are N\ — N 2 = 500, so N — Ni + 
N2 = 10 3 . With varying the size ratio oa/ci, the system 
volume V was changed such that the volume fraction of 
the soft-core regions defined by 



= (A^ + N 2 <rl)/V 



(2) 



was fixed at 0.9 mostly. We set <f> — 1 only in one case 
(in the lower panel of Fig. 8). With the mass ratio being 
rni/m% — (<7i/<72) 2 , we integrated the Newton equations 
using the leapfrog algorithm under the periodic boundary 
condition. The system temperature was controlled with 
the Nose-Hoover thermostat [H HI Ej. The time step 
of integration was 0.002r, where 



(7i y/mi/e. 



(3) 



Hereafter the time t. and the temperature T will be mea- 
sured in units of r and e/ks, respectively. 

We first equilibrated the system in a liquid state at 
T = 2 in a time interval of 10 3 and then quenched it 
to a lower final temperature with further equilibration 
in a period of t cq = 1.1 x 10 4 |23, 0. There was no 
appreciable time evolution in the pressure, the energy, etc 
in the time region t > 4 x 10 3 (see Fig. 8 as an example) 
|24j . The particles were well mixed and no indication of 
phase separation was observed in the final time region. 

In our study, the size ratio was in the range 1 < 
oi /c2 < 1-4. We saw no tendency of phase separa- 
tion. If <7i / go is too large, phase separation will be de- 
tected [23, |2g • We show typical particle configurations 
in Fig.l at the final simulation time t = 1.2 x 10 4 for (a) 
cr 2 / CTl = 1.1, (b) 1.2, (c) 1.225, and (d) 1.4 at (f> = 0.9. 
The system length V 1 ' 2 is (a) 35.03, (b) 36.81, (c) 37.27, 
and (d) 40.55 in units of <7\. They represent (a) a crystal 
state with point defects, (b) and (c) polycrystal states, 
and (d) a glass state. 



B. Sixfold orientation order 

In Fig.l, a large fraction of the particles are enclosed 
by six particles even at u\jai = 1.4. The particle config- 
urations are remote from other ordered structures such 
as the square structure Q. Therefore, we consider de- 
viations from the hexagonal order. The local crystalline 
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FIG. 1: Particle configurations for (a) ai/(J2 = 1.1, (b) 1.2, 
(c) 1.225, and (d) 1.4 at <j> = 0.9 and T = 0.2 taken at the 
final simulation time t = 1.2 x 10 4 . Smaller (larger) circles 
represent the smaller (larger) particles. For the visualization 
purpose the diameters of the circles in the snapshots are taken 
as Ag\ and A02 with A < 1 for the two species. 



order is represented by a sixfold orientation order vari- 
able [ll| • For each particle j we define 



Xj 



^ exp[6i6jk], 



(4) 



/cGbondcd 



where the summation is over the particles "bonded" to 
the particle j. In our case, the two particles j G a and 
k G (3 are bonded, if their distance rjk = \rj — rt| is 
shorter than R a p = 1.2h<j a p @- The upper cut-off R a p 
is slightly longer than the first peak position of the pair- 
correlation function g a p{r). The Ojk is the angle of the 
relative vector rj — with respect to the x axis. For a 
perfect triangular crystal of a one component system, the 
complex numbers Xj are an equal to 6 exp(6ia) with a 
being the common angle of one of the crystal axes with 
respect to the x axis. In the presence of disorder, the 
absolute values \Xj\ are significantly different from 6 for 
particles around defects. It is convenient to define a local 
crystalline angle ctj in the range < otj < 7r/3 by 



*j=X j /\Xj\=e e 



(5) 



In Fig. 2, we show the snapshots of the angles aj (j = 
1, • • • , 10 3 ) for the same particle configurations in Fig.l. 
The color map is illustrated in Fig. 3. We can clearly 
see point defects, grain boundaries, and glassy particle 
configurations. In (b) the grain boundaries are localized, 
while in (c) they are percolated. Recently, using a 2D 
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FIG. 2: Snapshots of the angles ctj in Eq.(5) for (a) ai/a-2 = 
1.1, (b) 1.2, (c) 1.225, and (d) 1.4 with the the color map in 
Fig.3. The data are common to those in Fig.l. Changeover 
from crystal to glass occurs with polycrystal as an intermedi- 
ate state. 



model of block copolymers, Vega et al. |27| numerically 
studied the grain boundary coarsening to obtain pictures 
of the orientation angles similar to our Fig. 2, though their 
system corresponds to one component particle systems. 
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FIG. 3: Color map for the angle aj in Eq.(5) for particle j 
at the center, around which the crystal order is perfect. The 
gray circle is the bonded region. The vector (red arrow) 
makes an angle of ay — 40° with respect to the horizontal 
axis. The color of particle j is then blue. 



Disorder variable 



We next introduce a new variable representing the de- 
cree of disorder. In terms of the difference — $j be- 
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FIG. 4: Disorder variable Dj in Eq.(6) for (a) a^/ax = 1.1, 
(a') 1.15, (b) 1.2, (c) 1.225, (c') 1.25, and (d) 1.4 with (j> = 0.9 
and t — 1.2 x 10 4 . The particle configurations in (a), (b), (c), 
and (d) are common to those in Figs.l and 2. Here the color 
changes in the order of rainbow. 



tween the bonded particle pairs, we define 

D i = E i^--**i 2 

Unbonded 

= 2 [1 — cos 6 (ay — at) 



(6) 



Unbonded 



for each particle j. This quantity is called the disorder 
variable. If the thermal vibrations are neglected, Dj van- 
ishes in single-component perfect crystals and is nonvan- 
ishing around defects. It takes large values of order unity 
almost everywhere in highly frustrated glass states. See 
the comment (iv) in the last section for appropriateness 
of this variable in glass and liquid. 

In Fig. 4, snapshots of Dj are shown for (a) 02 /ci = 
1.1, (a') 1.15, (b) 1.2, (c) 1.225, (c') 1.25, and (d) 1.4 at 
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FIG. 5: Structure factor for (a') o-z/oi = 1.15, (b) 1.2, and 
(b') 1.225 in Fig. 4. Bragg peaks can be seen in (a) and (b), 
while it resembles to that in liquid for (c) (see (c) in Fig. 7). 






FIG. 6: Disorder variable Dj in Eq.(6) in a polycrystal state 
at ffa/ffj = 1.225 for (bl) t = 6 x 10 3 and (b2) t = 8 x 10 3 , 
in the same run giving the panel (b) in Fig. 4 at t = 12 x 10 3 . 
Comparison of these three snapshots indicates very slow time 
evolution of the grain boundaries. See points (bl), (b2), and 
(b) in Fig. 7 also. 



the final states of the simulation runs at t = 1.2 x 10 4 . 
Those of (a), (b), (c), and (d) are taken from the same 
particle configurations as in the corresponding panels of 
Figs.l and 2. The color of the particles varies in the or- 
der of rainbow, being violet for Dj = and red for the 
maximum of Dj. In Fig. 4, the maximum of Dj is 1.71, 
2.64, 3.59, 4.36, 4.34, and 4.58 in (a)-(d) in this order. In 
crystals with &i ja\ close to unity, a small number of de- 
fects can be detected as bright points as in (a) and (a'). 
In polycrystals, defects are accumulated to form grain 
boundaries detectable as bright closed curves enclosing 
small crystalline regions, as in (b) and (c). With further 
increasing 02 /01 , defects are proliferated and a large frac- 
tion of the particles are depicted as bright points. In the 
largest size ratio in (d) , most of the particles are in disor- 
dered configurations. With varying 02/01, this crossover 
occurs in a narrow range around 1.2. 

In Fig. 5, the structure factor of the number density 
n(r) = ni(r) + n,2(r) is written for (a) 02/0-1 = 1.15, 
(b) 1.2, and (c) 1.225 to confirm the abruptness of this 
crystal-glass crossover. The structure factor in (a) ex- 
hibits Bragg peaks showing translational order, while 
that in (c) is similar to that in liquid but still retains the 
sixfold angular symmetry. For the intermediate case (b), 
the sixfold symmetry is evidently present and the trans- 



lational order is being lost. (See Fig. 9 below for structure 
factors in typical cases far from the transitions.) We note 
that similar structure factors were taken from a quasi 2D 
colloid suspension around the melting [l3j . 

The timescale of the particle configurations becomes 
exceedingly slow in glass states 0, 0, HI lU IH3- Also 
in polycrystal states, the motions of the grain bound- 
aries become slow with increasing 0-2/01, while the grain 
boundaries coarsen to disappear in one component sys- 
tems on a rapid timescale (see the corresponding curve 
in Fig. 7) [22]] • In Fig. 6, we present two additional snap- 
shots of Dj(t) at 02/01 = 1.225 for t = 6 x 10 3 and 
t = 8 x 10 3 , while the panel (b) in Fig. 4 is the snapshot 
at t = 12 x 10 3 in the same run. These three snapshots 
exhibit percolated grain boundaries with only small dif- 
ferences on large scales, indicating pinning of the grain 
boundaries. The panel (b) in Fig. 2 demonstrates that the 
system is a polycrystal. As a result, we cannot deduce 
the life time of the gran boundaries from our simulation 
in this case. 
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FIG. 7: Relaxation of the disorder parameter D(t) in Eq.(7) 
at <t> = 0.9 for 02/01 = 1.0, 1.175, 1.2, 1.225, and 1.4 from 
below. The temperature T is lowered from 2 to 0.2 at t — 10 3 . 
For the one component case (72/0-1 = 1.0, it takes place on a 
timescale of 50, as shown in the expanded inset. With size 
dispersity, D(t) relaxes slowly on timescales of order 4 x 10 3 . 
Snapshots of Dj(t) at two points (bl) and (b2) on the curve 
of 02/01 = 1.2 are given in Fig. 6. 



D. Degree of overall disorder 

We now introduce a single "disorder parameter" rep- 
resenting the degree of overall disorder by taking the av- 
erage over all the particles, 



D(t) 



where the time-dependence of Dj(t) and D(t) is explicitly 
written. In Fig. 7, we show time evolution of D{t), where 
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FIG. 8: Disorder parameter D in Eq.(8) for T = 0.2 and 2, 
with <j> — 0.9 in the upper panel and 1 in the lower panel. 
Liquid states are realized on the curve of T = 2 in the upper 
panel. The other curves show the crystal-glass crossover. 



quenching is from liquid at t = 10 3 . It undergoes very 
slow time evolution with finite size dispersity in polycrys- 
tal and glass, in accord with Fig. 6. For o^/fi = 1.4, D(t) 
increases upon quenching (see Fig. 9 below for its reason). 
In the time region t > 4x 10 3 , we can see no appreciable 
relaxation in these curves. For the one component case, 
the relaxation from liquid to crystal terminates rapidly 
on a timescale of 50. 

However, the curves in Fig. 7 with size dispersity 
weakly depend on time around the average even in ap- 
parent steady states. For example, D(t) is 1.10 in (bl) of 
Fig.6, 1.17 in (b2) of Fig.6, and 1.11 in (b) of Fig.4. This 
temporal fluctuations should diminish for larger system 
size. Its deviation from the time average became largest 
when the grain boundaries appreciably moved in poly- 
crystal states. For each simulation run, we defined the 
time average of D{t) as 



D 



1 

id 



dtD(t), 



(8) 



tf—td 



where if(= 1.2 x 10 4 ) is the terminal time of the sim- 
ulation run and id(= 2 x 10 3 ) is the width of the time 
interval of taking data. We regard D as a steady state 




FIG. 9: Disorder parameter D in Eq.(8) as a function of <Jija\ 
and T for <j> = 0.9. Typical structure factors for (A) crystal, 
(B) glass, and (C) liquid in the upper panel, where the corre- 
sponding points are indicated in the lower panel. 



average though glass states may further relax on longer 
timescales. 

In Fig. 8, we plot I) as a function of the size ratio. In 
the upper panel, where <f> = 0.9, liquid states are realized 
for any (J%I<J\ at T — 2, while the system is crystalline 
for a%Jcr\ < 1.2 and glassy for larger cr^/cri at T = 0.2. 
In the range o^/fi > 1.3, D in the liquid state at T = 2 
becomes smaller than D in the glass state at T = 0.2. 
This is because D increases weakly with increasing 02 /01 
in liquid and increases more strongly in glass. In the 
lower panel, where <f> — 1, the system crosses over from 
crystal to glass both for T — 0.2 and 2, and D increases 
rather abruptly around o^/ci ~ 1.2. For cr-z/ai < 1.1, D 
takes a small positive number due to the thermal motions 
of the particles. In Fig. 9, D is plotted as a function of 
02/ci and T. It shows the overall behavior of D. That 
is, D is small in crystal and increases abruptly in glass 
and liquid. Interestingly, for (J-xjuy > 1.25, D decreases 
with increasing T from glass to liquid (see Fig. 7). For 
such size ratios, highly disordered particle configurations 
can be pinned at low T and the thermal motions at high 
T can relax them. 



E. Defect-mediated melting 

In our simulations at fixed density, we observed defect 
proliferation at the melting (as well as at the crystal-glass 
crossover) and no coexistence of crystal and liquid regions 
separated by sharp interfaces. The system became highly 
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FIG. 10: Polycrystal configurations of otj in Eq.(5) (top) and 
Dj in Eq.(6) (bottom) at t = 6x 10 3 (left) and 6.2 x 10 3 (right) 
in the one component case az/ai = 1 at T = 1.3 and = 0.9. 
The system is intermediate between crystal and liquid. The 
defect structure is evolving rapidly on a timescale of 50. 



heterogeneous (as in Fig. 10 below), but no nucleation 
process could be detected. Figure 9 shows that D changes 
continuously along the T axis at each <Jij<7\ including 
the one component limit aija\ — 1. Similarly, in a 2D 
Lenard-Jones system with N — 256 at <j> = 0.8, Frenkel 
and Mctague detected no discontinuity in the average 
pressure and energy Theoretically, the 2D melting 
can be either continuous or first ordrer depending the 
specific details of the system [28|, |29J. It is a delicate 
problem to determine its precise nature in the presence 
of the heterogeneity developing at the transition [l3L FLU 

BG3. 

To visualize the physical process involved at the melt- 
ing, we display snapshots of otj (t) and Dj (t) at t = 6000 
and 6200 in Fig. 10 in the one component case at T = 1.3, 
where the change of D is abrupt in Fig. 9. The D(t) in 
Eq.(7) is 1.80 at t = 6000 and 1.60 at t = 6200. We 
can see percolated grain-boundary patterns and chains of 
point defects. The area fraction of the crystalline regions 
with small Dj continuously decreases (increases) with 
further raising (lowering) the temperature. We mention a 
simulation by McTague et al. in a one component system 
with soft-disk r -6 potentials 30], reporting the presence 
of both free dislocations and many grain boundaries at 
the melting. Some authors already pointed out relevance 
of grain boundaries in the 2D melting [2^. l3l| . Using 
inherent-structures theory, Somer et al. 32] found per- 
colated grain boundaries in "inherent structures" after a 



hexatic-to- liquid transition. Among many experiments, 
grain boundaries were evidently shown in Ref. 15]. 

We notice close similarity between the snapshots of the 
polycrystal states in Figs. 6 and 10. However, very differ- 
ent are the timescales of the dynamics of Dj(t) without 
and with size dispersity. Indeed, the patterns in Fig. 10 
changed appreciably on a rapid timescale of 50, while 
the large scale patterns in Fig. 6 were nearly frozen in our 
simulation time. 



III. SUMMARY AND REMARKS 

In summary, using MD simulation on a 2D LJ binary 
mixture, we have investigated the effects of the size dis- 
persity in the range 1 < <J\j 02 < 1.4 and the temperature 
in the particle configurations at fixed average density and 
composition. Our main objective has been to visualize 
defects, so the system size (N = 10 3 ) has been chosen 
to be rather small. Larger system sizes are needed to 
get reliable correlation functions of the density and the 
sixfold orientation variable. 

We summarize our main results and give remarks. 

(i) We have displayed the angle variable otj defined by 
Eq.(5) in Fig. 2 and the disorder variable Dj(t) defined 
by Eq.(6) in Fig.4 at low T = 0.2. The snapshots of 
these variables evidently show how the particle config- 
urations become disordered with increasing the size ra- 
tio. Those of Dj(t) provide the real space pictures of 
the defect structures on various spatial scales. We find 
polycrystal states with grain boundaries between crystal 
and glass. The motions of the grain boundaries are much 
slowed down with size dispersity, as in Figs. 6 and 7. 

(ii) The disorder parameter D(t) in Eq.(7) or its time av- 
erage D in Eq.(8) is a measure of overall disorder. As in 
Fig. 7, the relaxation of D(t) after quenching from a high 
to a low temperature occurs on a very long timescale with 
size dispersity, while it relaxes much faster in one com- 
ponent systems. The steady state average D is small in 
crystal and increases abruptly in glass and liquid, with 
increasing 01/02 or T, as in Figs. 8 and 9. 

(iii) In our system, the crystal-glass and crystal-liquid 
crossovers proceeded with increasing the defect density 
without nucleation. Remarkable resemblance is note- 
worthy between the polycrystal patterns of Dj (t) at the 
crystal-glass transition in Figs. 4 and 6 and those in the 
one component case at the crystal-liquid transition in 
Fig. 10. However, the timescale of the defect structure is 
drastically enlarged with increasing 01/02. In these two 
transitions, the disorder parameter D increases abruptly 
but continuously from small (crystal) values to large 
(glass or liquid) values. In these cases, polycrystal states 
appear with large scale heterogeneities in Dj(t), as can 
be seen in Figs. 4, 6, and 10. 

(iv) The particle configurations will increasingly deviate 
from the hexagonal order in the crossover from crystal 
to glass or liquid. They might become rather closer to 
other ordered structures for some fraction of the particles 
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HQ. In such cases, large values of Dj(t) and D will have 
only qualitative meaning, since they represent deviations 
from the hexagonal order. 

(v) We should study the pinning mechanism of grain 
boundaries in polycrystal in the presence of size disper- 
sity. We should also examine the dynamical properties 
such as the diffusion constant, the shear viscosity, and 
the time-correlation functions for various degrees of dis- 
order. They have been calculated around the liquid-glass 
transition |E 0, II ©, G3 

(vi) For the pair potentials in Eq.(l) and for our limited 
simulation time, we have detected no tendency of phase 
separation. By increasing the repulsion among the differ- 
ent components, we could study nucleation of crystalline 
domains in a glass matrix, for example. 

(vii) We will report shortly on the shear flow effect at the 



crystal-glass and crystal-liquid transitions in 2D. It has 
already been studied at the liquid-glass transition [SillCj. 
It is of interest how an applied shear affects the defect 
structure and induces plastic deformations |33| . 
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